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Abstract 

This paper addresses the question of how to modify in aerodynamic design to 
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putational feasibility of using control theory for such a purpose. 
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This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18107 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 


i 


1. Introduction and historical survey 


Computers have had a twofold impact on the science of aerodynamics. On the one hand 
numerical simulation may be used to gain new insights into the physics of complex flows. On the 
other hand computational methods can be used by engineers to predict the aerodynamic 
characteristics of alternative designs. Assuming that one has the ability to predict the 
performance, the question then arises of how to modify the design to improve the performance. 

This paper is addressed to that question. 

Prior to 1960 computational methods were hardly used in aerodynamic analysis. The 
primary tool for the development of aerodynamic configurations was the wind tunnel. Shapes 
were tested and modifications selected in the light of pressure and force measurements together 
with flow visualization techniques. Computational methods are now quite widely accepted in the 
aircraft industry. This has been brought about by a combination of radical improvements in 
numerical algorithms and continuing advances in both speed and memory of computers. 

If a computational method is to be useful in the design process, it must be based on a 
mathematical model which provides an appropriate representation of the significant features of the 
flow, such as shock waves, vortices and boundary layers. The method must also be robust, not 
liable to fail when parameters are varied, and it must be able to treat useful configurations, 
ultimately the complete aircraft. Finally reasonable accuracy should be attainable at reasonable 
cost. Much progress has been made in these directions [1—10]. In many applications where the 
flow is unseparated, including designs for transonic flow with weak shock waves, useful predictions 
can be made quite inexpensively using the potential flow equation [1—4]. Methods are also 
available for solving the Euler equations for two— and three-dimensional configurations up to a 
complete aircraft [5-10]. Viscous simulations are generally complicated by the need to allow for 
turbulence: while the Reynolds averaged equations can be solved by current methods, the results 
depend heavily on the choice of turbulence models. 
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Given the range of well proven methods now available, one can distinguish objectives for 
computational aerodynamics at several levels: 

1) Capability to predict the flow past an airplane or important components in different flight 
regimes such as take-off or cruise, and off design conditions such as flutter. 

2) Interactive design calculations to allow rapid improvement of the design. 

3) Automatic design optimization. 

Substantial progress has been made toward the first objective, and in relatively simple 
cases such as an airfoil or wing in inviscid flow, calculations can be performed fast enough that the 
second objective is within reach. The third objective has also been addressed for various special 
cases. In particular it has been recognized that the designer generally has an idea of the kind of 
pressure distribution that will lead to the desired performance. Thus it is useful to consider the 
problem of calculating the shape that will lead to a given pressure distribution. Such a shape does 
not necessarily exist, unless the pressure distribution satisfies certain constraints, and the problem 
must therefore be very carefully formulated: no shape exists, for example, for which stagnation 
pressure is attained over the entire surface. 

The problem of designing a two dimensional profile to attain a desired pressure distribution 
was first studied by Lighthill, who solved it for the case of incompressible flow by conformally 
mapping the profile to a unit circle [11]. The speed over the profile is 


q = <f> e l h (1.1) 

where <p is the potential for flow past a circle, and h is the modulus of the mapping function. The 
solution for <f> is known for incompressible flow. Let qd be the desired surface speed. Then the 
surface value of h can be obtained by setting q = qd in equation (1.1), and since the mapping 
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function is analytic, it is uniquely determined by the value of h on the boundary. A solution 
exists for a given speed at infinity only if 


1 

Si 


(b qd0 = 


(1.2) 


and there are additional constraints on q if the profile is required to be closed. 


Lighthill’s method was extended to compressible flow by McFadden [12]. Starting with a 
given shape, and a corresponding mapping function h^, the flow equations can be solved for the 
potential <p^\ which now depends on h^. A new mapping function h^ is then determined by 
setting q = qd in equation (1.1), and the process is repeated. In the limiting case of zero Mach 
number the method reduces to Lighthill’s method, and McFadden gives a proof that the iterations 
will converge for small Mach numbers. He also extends the method to treat transonic flow 
through the introduction of artificial viscosity to suppress the appearance of shock waves, which 
would cause the updated mapping function to be discontinuous. This difficulty can also be 
overcome by smoothing the changes in the mapping function. Such an approach is used in a 
computer program written by the author for Grumman Aerospace. It allows the recovery of 
smooth profiles that generate flows containing shock waves, and it has been used to design 
improved blade sections for propellers [13]. A related method for three dimensional design was 
devised by Garabedian and McFadden [14]. In their scheme the steady potential flow solution is 
obtained by solving an artificial time dependent equation, and the surface is treated as a free 
boundary. This is shifted according to an auxiliary time dependent equation in such a way that 
the flow evolves toward the specified pressure distribution. 


Another way to formulate the problem of designing a profile for a given pressure 
distribution is to integrate the corresponding surface speed to obtain the surface potential. The 
potential flow equation is then solved with a Dirichlet boundary condition, and a shape correction 
is determined from the calculated normal velocity through the surface. This approach was first 
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tried by Tranen [15]. Volpe and Melnik have shown how to allow for the constraints that must be 
satisfied by the pressure distribution if a solution is to exist [16]. The same idea has been used by 
Henne for three-dimensional design calculations [17]. 

The hodograph transformation offers an alternative approach to the design of airfoils in 
transonic flows. Garabedian and Korn achieved a striking success in the design of airfoils to 
produce shock— free transonic flows by using the method of complex characteristics to solve the 
equations in the hodograph plane [18]. Another design procedure has been proposed by Giles, 
Drela and Thompkins [19], who write the two-dimensional Euler equations for inviscid flow in a 
streamline coordinate system, and use a Newton iteration. An option is then provided to treat the 
surface coordinates as unknowns, while the pressure is fixed. 

Finally, Hicks and Henne have explored the possibility of meeting desired design objectives 
by using constrained optimization [20]. The configuration is specified by a set of parameters, and 
any suitable computer program for flow analysis is used to evaluate the aerodynamic 
characteristics. The optimization method then selects values of these parameters that maximize 
some criterion of merit, such as the lift— to— drag ratio, subject to other constraints such as 
required wing thickness and volume. In principle this method allows the designer to specify any 
reasonable design objectives. The method becomes extremely expensive, however, as the number 
of parameters is increased, and its successful application in practice depends heavily on the choice 
of a parametric representation of the configuration. 

The purpose of this paper is to propose that there are benefits in regarding the design 
problem as a control problem in which the control is the shape of the boundary. A variety of 
alternative formulations of the design problem can then be treated systematically by using the 
mathematical theory for control of systems governed by partial differential equations [21]. 

Suppose that the boundary is defined by a function f(x), where x is the position vector. As in the 
case of optimization theory applied to the design problem, the desired objective is specified by a 
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cost function I, which may, for example, measure the deviation from a desired surface pressure 
distribution, but could also represent other measures of performance such as lift and drag. The 
introduction of a cost function has the advantage that if the objective is unattainable, it is still 
possible to find a minimum of the cost function. Now a variation in the control 6f leads to a 
variation £1 in the cost. It is shown in the following sections that SI can be expressed to first order 
as an inner product of a gradient function g with 81: 

+ 

<51 = (g , Si) 

Here g is independent of the particular variation 61 in the control, and can be determined by 
solving an adjoint equation. Now choose 

<5f = -Ag 

where A is a sufficiently small positive number. Then 

SI = -A(g,g) < 0 

assuring a reduction in I. After making such a modification, the gradient can be recalculated and 
the process repeated to follow a path of steepest descent until a minimum is reached. In order to 
avoid violating constraints, such as a minimum acceptable wing thickness, the steps can be taken 
along the projection of the gradient into the allowable subspace of the control function. In this 
way one can devise design procedures which must necessarily converge at least to a local 
minimum, and which might be accelerated by the use of more sophisticated descent methods. 
While there is a possibility of more than one local minimum, the cost function can be chosen to 
reduce the likelihood of difficulties caused by such a contingency, and in any case the method will 
lead to an improvement over the initial design. The mathematical development resembles in 
many respects the method of calculating transonic potential flow proposed by Bristeau, 
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Pironneau, Glowinski, Periaux, Perrier and Poirier, who reformulated the solution of the flow 
equations as a least squares problem in control theory [4]. 

In order to illustrate the application of control theory to design problems in more detail, 
the following sections present design procedures for three examples. Section 2 discusses the design 
of two dimensional profiles for compressible potential flow when the profile is generated by 
conformal mapping. This leads to a generalization of the methods of Lighthill and McFadden. 
Section 3 discusses the same problem when the flow is governed by the inviscid Euler equations. 
Finally, Section 4 addresses the three dimensional design problem for a wing, assuming the flow to 
be governed by the inviscid Euler equations. The procedures which are presented require the 
solution of several partial differential equations at each step. The question of the most efficient 
discretization of these equations is deferred for future investigation. 
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2. Design for potential flow using conformal mappin; 

Consider the case of two dimensional compressible inviscid flow. In the absence of shock 
waves an initially irrotational flow will remain irrotational, and we can assume that the velocity 
vector <1 is the gradient of a potential <p. In the presence of weak shock waves this remains a fairly 
good approximation. Let £ , T and S denote vorticity, temperature and entropy. Then according 
to Crocco’s Theorem, vorticity in steady flow is associated with entropy production through the 
relation 

gx{+T?S=0 

Thus, the introduction of a potential is consistent with the assumption of isentropic flow, and 
shock waves are modelled by isentropic jumps. Let p, p, c and M be the pressure, density, speed 
of sound and Mach number q/c. Then the potential flow equation is 

V-pV0= 0 (2.1) 

where the density is given by 

P = { 1 + # M 2 (m 2 )} 1 ^ 1 (2.2) 

while 

c2 = * (2 ' 3) 

OD 

Here M is the Mach number in the free stream, and the units have been chosen so that p and q 
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have the value unity in the far field. Equation (2.2) is a consequence of the energy equation in the 
form 



Figure 1 


Suppose that the domain D exterior to the profile C in the z plane is conformally mapped 
onto the domain exterior to a unit circle in the a plane as sketched in Figure 1. Let R and 6 be 
polar coordinates in the a plane, and let r be the inverted radial coordinate 1/R. Also let h be the 
modulus of the derivative of the mapping function 


h = 


dz 


Now the potential flow equation becomes 

% { p ^ q ) + r ( r P^ r ) = 0 * n D 


(2.4) 


where the density is given by equation (2.1), and the circumferential and radial velocity 
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components are 


while 


J't’e 
u IT ’ 



(2.5) 


q 2 = u 2 + v 2 (2.6) 

The condition of flow tangency leads to the Neumann boundary condition 

v = E77r’ = ^ on ^ (2-^) 

In the far field the potential is given by an asymptotic estimate, leading to a Dirichlet boundary 
condition at r = 0 [2]. 

Suppose that it is desired to achieve a specified velocity distribution q^ on C. Introduce 
the cost function 


I = H (q-<5 d ) 2 it (2 8) 

c 

The design problem is now treated as a control problem where the control function is the mapping 
modulus h, which is to be chosen to minimize I subject to the constraints defined by the flow 
equations (2.2 — 2.7). 

A modification 6h to the mapping modulus will result in variations 6(f), £u, 6v and 6p to the 
potential, velocity components and density. The resulting variation in the cost will be 


<51 = J (q— q d ) «q it 


(2.9) 
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where on C q = u. Also 


<5u = r 


6< t>9 

TT" U 


6h 

T 


6v = r 2 


s<f> T 

1 T“ v 


1 


while according to equations (2.2) and (2.6) 



Hence 


6p = — ^ ( u< ^ u + vtfv) 

= p \ IT “ H ^ nS(f> 9 + vr ^ r ) 
c c 


It follows that 6(p satisfies 


L6<p - - ^ {pWp<t>Q ^jj) - r ^ (pM 2 r0 r -^) 


where 


L- d 
L = W 



d 

W 


pnv 


r l} 



d pnv d } 


( 2 . 10 ) 


Then if ip is any periodic differentiable function which vanishes in the far field 



pM 2 V0 • V^dS 


( 2 . 11 ) 


where dS is the area element rdrd0, and the right hand side has been integrated by parts. 


11 


Now we can augment equation (2.9) by subtracting the constraint (2.11). The auxiliary 
function ip then plays the role of a Lagrange multiplier. Substituting for <5q and integrating the 
term 


f (q-q d ) r ^0d0 

c 


by parts, we obtain 


- f (q — q d ) q ^ d0 


J 




q - q d 


1 d$- f ^L6<pdS + f 

D r D 


pM z V<p 


V^^dS 


Now suppose that ip satisfies the adjoint equation 


Lip = 0 in D 


( 2 . 12 ) 


with the boundary condition 



r q ~ q <n 


on C 


(2.13) 


Then integrating by parts 



L 6<p dS = - 


f pip T 6<pd9 
C 


(q - q d ) q + 


pM 2 v0 • v^ds 
t> 


and 


<51 = - 


(2.14) 
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Here the first term represents the direct effect of the change in the metric, while the area integral 
represents a correction for the effect of compressibility. 

Equation (2.14) can be further simplified to represent £1 purely as a boundary integral 
because the mapping function is fully determined by the value of its modulus on the boundary. 
Set 


l°g^ = f+i/3 


where 


f = log 


dz 

3 ? 


= log h 


and 


81 = 


<5h 

T 


Then f satisfies Laplace’s equation 


Af = 0 in D 


and if there is no stretching in the far field, f -* 0. Thus 

A Si = 0 in D 


and 81 -* 0 in the far field. 


Introduce another auxiliary function P which satisfies 


AP = *>M 2 V^ • in D 


( 2 . 15 ) 
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and 


P = 0 on C 


Then the area integral in equation (2.14) is 


J AP<5f dS = J 
D C 


^l? d * - | PA5fdS 


(2.16) 


and finally 


where 



6 = ^-(Q-q d )<l 


(2.17) 


(2.18) 


This suggests setting 


<5f=-Ag 

so that if A is a sufficiently small positive number 

<51 = -A J g 2 d0 < 0 

Arbitrary variations cannot, however, be admitted. The condition that f -» 0 in the far 
field, and also the requirement that the profile should be closed, imply constraints which must be 
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satisfied by f on the boundary C. Suppose that log 


dz' 


is expanded as a power series 


log 


dz' 


- I -S 

n=0 a 


(2.19) 


where only negative powers are retained because otherwise ^ would become unbounded for large 
a. The condition that f -* 0 as a -* od implies 


c 0 - ° 


Also the change in z on integration around a circuit is 


Az = <J> 


^§ d <7 = 2iriCi 


so the profile will be closed only if 


c 1 = 0. 


On C equation (2.19) reduces to 


*C + * = \ ( a n cos + t> n sin n0) + i ^ (b Q cos n9 - a n sin n 6) 

n=0 n=0 


Thus a_ and b are the Fourier coefficients of L-,, and these constraints reduce to 
n h 


a Q = 0, a 1 = 0, bj = 0 


In order to satisfy these constraints we can project g on to the admissible subspace for f n 
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by setting 


g = g - Aq - cos 6 - B 1 sin 6 (2.20) 

where 

A 0 = 3? j gd6 
C 

A 1 = ± J geos 0d0 (2.21) 

C 

B i = M 6 sin ede 

C 

Then 

‘(g-g) gd0 = 0 

C 


and if we take 


ff=-Ag 

it follows that to first order 


<51 -A 


g gd* = -A 


(g + g-g) gd*=-A 



If the flow is subsonic this procedure should converge toward the desired speed distribution 
since the solution will remain smooth, and no unbounded derivatives will appear. If, however, the 
flow is transonic, one must allow for the appearance of shock waves in the trial solutions, even if 
q^ is smooth. Then q — is not differentiable. This difficulty can be circumvented by a more 
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sophisticated choice of the cost function. Consider the choice 


-ij 


A, S z + A r 


dS 




dd 


( 2 . 22 ) 


where A^ and A 2 are parameters, and the periodic function S(0) satisfies the equation 


AjS- 


d 2s _ 

dr d 


(2.23) 


Then 


SI 


(Aj S <5S + ^2 SS)dO 


= | S (A l( 5S-A 2 ^£S)d0 


S<5qd0 


Thus S replaces q — q^ in the previous formulas, and if one modifies the boundary condition (2.13) 
to 



fSl 

E 


on C 


(2.24) 


the formula for the gradient becomes 


g = 



-Sq 


(2.25) 
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instead of equation (2.18). Then one modifies f by a step — Ag in the direction of the projected 
gradient as before. 


The final design procedure is thus as follows. Choose an initial profile and corresponding 
mapping function f. Then 

1) Solve the flow equations (2.2 — 2.7) for <p, u, v, q, p. 

2) Solve the ordinary differential equation (2.23) for S. 

3) Solve the adjoint equation (2.12) for ip subject to the boundary condition (2.24). 

4) Solve the auxiliary Poisson equation (2.15) for P. 

5) Evaluate 


„ _ d? Qn 
g-^-Sq 


on C, and find its projection g onto the admissible subspace of variations according to 
equations (2.20) and (2.21). 

6) Correct the boundary mapping function f^ by 



and return to step 1. 


3. Design for the Euler equations using conformal mapping 


This section treats the case of two dimensional compressible flow where the potential flow 
equation is replaced as a mathematical model by the inviscid Euler equations. Let p, p , u, v, E 
and H denote the pressure, density, Cartesian velocity components, total energy and total 
enthalpy. For a perfect gas 

P = (r-l) />{E -j(u 2 + v 2 )} (3.1) 

and 

pH = pE + p (3.2) 

where 7 is the ratio of specific heats. The Euler equations may then be written as 

< 3 - 3 > 

where x and y are Cartesian coordinates, t is the time coordinate and 



p 


'pu 

n 


pv 

w = 

pu 

, * = 

pu + p 

. g = 

pvu 

0 


pv 


puv 


pv + p 


pE 


puH 


pv H 


As in the previous section, suppose that the domain D exterior to the profile C in the z 
plane is mapped conformally onto the domain exterior to a unit circle in the a plane (see Figure 
1). Assume also that the outer boundary B of the domain is very far from the profile. Let the 
derivative of the mapping function be 
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Also let r and 9 be polar coordinates in the a plane, where in this case it is more convenient to 
take r as the true radial coordinate denoted by R in the previous section, and 9 is measured in the 
clockwise direction. Define the rotation parameters 


c = cos (0 - 9) , 


s = sin (0 — 9) 


(36) 


and rotated velocity components 



(3.7) 


Then the Euler equations become 

It (rh 2 w) + (hF) + (rhG) = 0 


(3.8) 


where 



’ p 


- P u 


pV 

w = 

pu 

py 

pE 

F = 

pU u + sp 
pU v — cp 
pUH 

G = 

pVu + cp 
pVv + sp 
.PVH 


(3.9) 


Then the flow is determined as the steady state solution of equations (3.8) and (3.9), subject to 
the flow tangency condition 


V = 0 onC 


(3.10) 


At the far field boundary B conditions can be specified for incoming waves, while outgoing waves 
are determined by the solution. 
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In contrast to the case of potential flow, the pressure is not determined solely by the speed, 
and assuming that one wishes to control the surface pressure distribution, a suitable cost function 
is 




(P-P d )' 


it 


(3.11) 


where p d is the desired pressure. A modification to the mapping function will influence equations 
(3.8) and (3.9) through changes $i and 60 in both the modulus and argument of ^ , finally 
leading to a variation in the cost function 


a = f (p - P d ) 6p dO 

C 


(3.12) 


where 6p is the variation of the pressure. 


Now the mapping variations cause variations in the rotation parameters 


6s = c 60 , 

6( hs) = s6h 4- he 60 , 


6c = - s 60 
6{ he) = c6h — hs 60 


(3.13) 


Define the Jacobian matrices 


A = 


di 


B = fe 

<7W 


C = sA - cB , 


D = cA + sB 


(3.14) 
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Then the variation 6w in w satisfies 


!j(hCtfw) + ^(rhD«w) = -|j(F & + hG«/3) -|j r(G<h - hF$) 


(3.15) 


Also 


5V = 0 on C 


(3.16) 


At the outer boundary there will be no variation in characteristic variables corresponding to 
incoming waves. If we take the outer boundary B at a fixed radius, incoming waves correspond to 
negative eigenvalues of D. Suppose that D is represented as T A T — \ where A is a diagonal 
matrix containing its eigenvalues, and the columns of T are eigenvectors of D. Define 

5w = T — 6w 


and <$w as the components of <5w corresponding to negative eigenvalues of T. Then 


6w =0 on B 


(3.17) 


Since 6w satisfies the constraint (3.15), we can replace equation (3.12) by 


51 = 


(p - Pd)£pd0 — (hC<5w) + ^ (rhDtfw) jdrd# 

D 


f {Is ( Fa + M3$) + If r(Ga - hF«/5)}drd0 


(3.18) 


where the vector ^ is a Lagrange multiplier, and the superscript T denotes the transpose. Suppose 
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that ip is the steady state solution of the adjoint equation 

$-C T f£-rD T $=0 inD (3.19) 

At the outer boundary B conditions can be specified for incoming waves, corresponding to positive 
T _i ^ m 

eigenvalues of D = T AT. Define 


$ = T T ^- 


and as the components of $ corresponding to positive eigenvalues of D. Then we can set 


ijj+ = 0 on B 


(3.20) 


If we integrate equation (3.18) by parts the contribution 


rh^> D£wd# = rh$A£wd0 




vanishes because of the complementary boundary conditions (3.17) and (3.20) satisfied by £w and 
xp at the outer boundary. If <$h and 6(3 decay fast enough in the far field the contribution 

|r(G£h-hF£/?) d0 
will also be negligible. Thus we find that 


<51 = J (p — P d ) 6pd0 


+ 


i> L (hDtfw + G6h - hF 60) d0 + J 
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where 

J = J J {(F T ^ + G T rtf- r )fh + (G T ^ - F T r^ r ) h£/?}drd0 (3.21) 

D 


Also 


hD5w + GSh - hFSp = 6(hG) = h 


'0 ' 


'0 

c6p 

+ P 

ff(hc) 

si5p 

6( hs) 

. 0 


_ 0 


Thus using the relations (3.13) 


61 = 


(p-p^tfpd# + (c^ + s^^phd# + f p(cV»2 + sip^)6hd0 - p(s^ 2 ~ c , ip^)6phdd + J 

C C C 


Now let ip satisfy the boundary condition 


h(c^ 2 + s ip 3 ) = - (p - p d ) on C 


(3.22) 


Then 



(P”P d ) P T^ - 

6 


(s^ _ c^g)ph<5/M0 + 


J 


(3.23) 


Finally we can use the fact that the mapping function is fully determined by its boundary 
value to reduce J to a boundary integral. Set 

>°sa §= f+i » 


where 


f = log 


dz 


= logh 
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( 3 . 24 ) 


( 3 . 25 ) 


( 3 . 26 ) 


( 3 . 27 ) 
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Then 


J = J J (AP & + AQ 60 ) drd# 
D 



& jfr + hp(s^2 + ^3) — Q ^ W d 0 


— | 7J7 + ^P( s ^2 + ~^W 


Thus finally 


SI = 



(3.28) 


where 

SlT> OQ 

6 = W + W - ( p “ p d^ p 


(3.29) 


As in the previous section, an appropriate modification of f is 


Sf = — Ag 


where g is the projection of g onto the admissible subspace of variations defined by equations 
(2.20) and (2.21), and A is a sufficiently small positive number. Then 

SI = -A [ g 2 d0 < 0 
C 


If the flow is transonic, shock waves are likely to be formed, and again it may be desirable 
to use a more sophisticated cost function to produce a smooth shape change. In this case we can 
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set 


I = 



AjS 2 + A 2 


dS] 

Kd 


dO 


(3.30) 


where Ai and A 2 are positive parameters, and the periodic function S(0) satisfies the equation 


A 



d 2 S 




(3.31) 


Then 


a = f (A lS «s + A 2 ^^ffi)d9 

c 


J S^SS - 6S) d8 


= [ S<5pd0 
C 


Thus S replaces p — p^ in the previous formulas. If one modifies the boundary condition (3.22) to 


h(c ^2 + s^g) = - S on C 


(3.32) 


the formula for the gradient becomes 


Z~W + 


-Sp 


(3.33) 


instead of equation (3.29), and an appropriate modification of f is again — Ag. 
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The final design procedure using the Euler equations is thus as follows. Choose an initial 
profile and corresponding mapping function f. Then 

1) Solve the flow equation (3.8) for w by integrating to a steady state. 

2) Solve the ordinary differential equation (3.31) for S. 

3) Solve the adjoint equation (3.19) with the boundary conditions (3.20) and (3.32) for ip by 
integrating to a steady state. 

4) Solve the auxiliary Poisson equations (3.24) and (3.26) for P and Q. 

5) Evaluate 



on C, and find its projection g onto the admissible subspace of variations according to 
equations (2.20) and (2.21). 

6) Correct the boundary mapping function f^ by 

«=- Ag 


where A > 0, and return to step 1. 
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4. Wing design using the Euler equations 


In order to illustrate further the application of control theory to aerodynamic design 
problems, this section treats the case of three-dimensional wing design, again using the inviscid 
Euler equations as the mathematical model for compressible flow. In this case it proves 
convenient to denote the Cartesian coordinates and velocity components by Xp Xg and Up u 
Ug, and to use the convention that summation over i = 1 to 3 is implied by a repeated index i. 
The three-dimensional Euler equations may then be written as 
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) 


dw 

cTT + 


5f i 

W. 

i 


= 0 


(4.1) 


where 
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(4.2b) 


Also 


P = (7“1)p(E 



> 


pH = pE + p 


(4.3) 


29 


Consider a transformation to coordinates Xp X2> X^ where 
dx. 


H ; • = 


ij cCC 


J = det(H), 


-i 


(44) 


The Euler equations can now be written as 


dW dF i _ o 
3T + 3X7" 0 


(4.5) 


where 


dX { 

W = Jw, F i = J^if 

J J 


(4.6) 


Define the contravariant velocity vector 
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U 2 
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(4.7) 


Then 


F i = J 


pUi 

axj 

^ u i u i + ^5q; p 
axj 

^ U i u 2 + ^ P 

^ U i u 3 + ^ P 
LpU 


(4.8) 


Assume now that the new coordinate system conforms to the wing in such a way that the wing 
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surface B w is represented by = 0. Then the flow is determined as the steady state solution of 
equation (4.5) subject to the flow tangency condition 


U 0 = 0 on B 
2 w 


(4.9) 


At the far field boundary, conditions can be specified for incoming waves as in the 
two-dimensional case, while outgoing waves are determined by the solution. 


Suppose now that it is desired to control the surface pressure by varying the wing shape. It 
is convenient to retain a fixed computational domain. Variations in the shape then result in 
corresponding variations in the mapping derivatives defined by H. Introduce the cost function 


i = z}J(p-P d ) 2 dx 1 dx, 

B w 


(4.10) 


where pd is the desired pressure. A variation in the shape will cause a variation £p in the pressure 
and consequently a variation in the cost function 


£ = JJ (P - P d ) 6p dX x dX 3 
B 


(4.11) 


w 


Since p depends on w through the equation of state (4.3), the variation <5p can be 
determined from the variation 6w. Define the Jacobian matrices 


A. = 


9 h 

0W 


°i “ Vi 


(4.12) 


Then 


jUj-tiFj )- 0 


(4.13) 
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where 


<§F. = C-<5w + 


axj 

-gr) fj 


(4.14) 


and for any differentiable vector ip 


dv 

b 1 


| n-^^Fj ds 

bouda r ies 


(4.15) 


where n^, and n^ are the components of a unit vector normal to the boundary. On the wing 
surface B w , n^ = n^ = 0 and it follows from equation (4.9) that 


' 0 


'0 

dX 0 

nq* 


dX 0 

dX 0 

V p 

+ P 

dX 0 

dX 0 

^ Sp 


dX 0 

. 0 


.0 


Suppose now that ip is the steady state solution of the adjoint equation 


(4.16) 



(4.17) 


At the outer boundary incoming characteristics for ip correspond to outgoing characteristics for 
Sw. Consequently, as in the two-dimensional case, one can choose boundary conditions for ip such 
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that 


n i^C'i 6w = 0 


If the coordinate transformation is such that <5(JH is negligible in the far field, the only 
remaining boundary term is 



ip T 6F 2 dX x dX 3 


Let ip satisfy the boundary condition 



onB w 


(4.18) 


Then, since it follows from equation (4.17) that 



6w dV = 0 


we find that 




+ iM 



p dXjdXg 


(4.19) 
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Figure 2 

A convenient way to treat a wing is to introduce sheared parabolic coordinates through the 
transformation 


x = ^ {x 2 - (Y + S(X,Z)) 2 } 
y = x(y + S(X,Z)) 


Z = Z 

» 

t 

Here x, y, z are Cartesian coordinates, and X and Y-fS correspond to parabolic coordinates 
generated by the mapping 

x + iy = \ (X + i(Y + S)) 2 

at a fixed span station Z. The surface Y=0 is a shallow bump corresponding to the wing surface, 
1 with a height S(X,Z) determined by the equation 


X + iS = 


2 (*s + '*,) 


where x (z) and y (z) are coordinates of points lying in the wing surface. We now treat S(X,Z) as 

S' S' 
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the control. 


In this case 


H = 


'X-(Y+S)S X 
Y + S + XS X 


0 


-(Y+S) -(Y+S)S Z 
X XS Z 

0 1 


while 

J = X 2 + (Y+S) 2 


and 
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(X— (Y+S)S X -JS Z 
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Also 

63 = 2(Y+S) £S 
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(<5SS X +(Y+S)£S X ) - (&J S Z +J^S Z ) 

0 
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Inserting these formulas in equation (4.19) we find that the volume integral in 61 is 

J ^ tfSfjdV 

D 

- £ {(ffi+xas^ + (6SS x +(Y+S)6S x )f 2 + (6JS z +J6S z f 3 )} dV 

+ | Vz &MV 

D 

where S and 6S are independent of Y. Therefore, integrating over Y, the variation of the cost 
function can be reduced to a surface integral of the form 

SI = (P(X,Z)6S + Q(X,Z)6S X + R(X,Z)6S z )dXdZ 


Also the shape change will be confined to a bounded region of the X-Z plane, so we can integrate 
by parts to obtain 


)I 


61= (P - 


- 1|) 6S dX dZ 


Thus to reduce I we can choose 


5Rx 


SS = -A(P-^--^) 


where A is sufficiently small and non— negative. 

In order to impose a thickness constraint we can define a baseline surface Sq(X,Z) below 
which S(X,Z) is not allowed to fall. Now if we take A = A(X,Z) as a non-negative function such 
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that 


S(X,Z) + <5S(X,Z) > S Q (X,Z) 

Then the constraint is satisfied, while 

SI = -JJ A(P-|^-||) 2 dXdZ < 0 
B 
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5. Conclusion 


The purpose of the last three sections is to demonstrate by representative examples that 
control theory can be used to formulate computationally feasible procedures for aerodynamic 
design. The cost of each iteration is of the same order as two flow solutions, since the adjoint 
equation is of comparable complexity to the flow equation, and the remaining auxiliary equations 
could be solved quite inexpensively. Provided, therefore, that one can afford the cost of a 
moderate number of flow solutions, procedures of this type can be used to derive improved 
designs. The approach is quite general, not limited to particular choices of the coordinate 
transformation or cost function, which might in fact contain measures of other criteria of 
performance such as lift and drag. For the sake of simplicity certain complicating factors, such as 
the need to include a special term in the mapping function to generate a corner at the trailing 
edge, have been suppressed from the present analysis. Also it remains to explore the numerical 
implementation of the design procedures proposed in this paper. 
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